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1. Introduction 



Modern accelerator experiments require precise perturbative calculations for the event rate 
of a variety of physical processes. Jets, electroweak gauge bosons and heavy quarks are 
being produced copiously at the Tevatron and the LHC. The precision of the measurements 
of physical masses, coupling parameters and the structure of colliding hadrons depends 
significantly on theoretical uncertainties which are better controlled at higher orders in 
perturbation theory The exclusion of hypotheses for novel particles and interactions is 
more significant when candidate signal processes are predicted accurately. With the arrival 
of new discoveries, the nature of physics laws will be deciphered more confidently with the 
aid of solid quantitative theory predictions. 

Our abilities to simulate complicated physical processes beyond the leading order (LO) 
have been improved dramatically in the last few years. At next-to-leading-order (NLO), 
previously inaccessible calculations with up to five particles in the final state are now 
possible [1]. Basic collider processes with fewer particles have also been computed at next- 
to- next-to-leading (NNLO) order in QCD [2-13]. For hadron colliders, the experimental 
frontier in particle physics, only cross-sections for 2 — > 1 processes have been computed 
at NNLO. Such computations must be extended to 2 — > 2 processes which are relevant to 
the experimental program. These include top and bottom quark pair production, inclu- 
sive jet production, electroweak diboson production, electroweak gauge boson and Higgs 
production in association with jets, single top production and beyond the Standard Model 
signals. It is unclear whether existing methods are suited to this task, and refinements 
of traditional methods or the development of new ones are required in order to face the 
increased complexity of such calculations. 

A fundamental technical difficulty in NNLO calculations is the appearance of multi- 
dimensional integrals over the momenta of up to two additional real or virtual particles 
with respect to the Born process. These integrals are separately infrared divergent and only 
their sum is finite. Higher order computations are performed almost exclusively within di- 
mensional regularization, where real and virtual corrections are expanded in a dimensional 
regulator e = 2 — ^ , where D is the number of dimensions. Laurent expansions in e are in- 
tricate in the presence of overlapping singularities. In this paper we present a new method 
for the calculation of the Laurent series in e of multidimensional integrals which typically 
appear in NNLO computations and generic higher order computations. 

Existing methods which tackle or bypass the problem of overlapping singularities are 
based on differential equations [14-16], Mellin-Barnes representations [17-20] and sector de- 
composition [21-23]. The first two approaches can be applied to the calculation of virtual or 
inclusive real radiation corrections. A subtraction method can reduce the problem of fully 
differential real radiation calculations at NNLO to inclusive phase-space calculations [10-12] 
enabling the method of differential equations and Mellin-Barnes for fully differential calcu- 
lations. Sector decomposition can be applied universally, for virtual inclusive phase-space 
integrations and fully differential integrations of real radiation matrix-elements. 

Sector decomposition has been the first successful method for performing fully differ- 
ential NNLO calculations for hadron collider processes [7,8]. This is largely attributed to 
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the conceptual simplicity of the method and its algorithmic nature which permits a full au- 
tomatization. The algorithm eliminates overlapping singularities by slicing the integration 
volume such that variables which contribute to an overlapping singular limit are ordered. 
In this way, the singularity is always factorized and it appears as a single singular limit 
of only the largest variable. While this algorithm leads to numerically stable evaluations, 
it requires a large number of integrals due to the slicing of the integration volume. This 
hinders the application of the method to processes with more complicated matrix-elements. 

We present here an alternative method for the factorization of overlapping singularities. 
We have observed that a factorization is possible by means of simple rescaling of the 
integration variables and non-linear transformations which preserve the geometry of the 
integration boundaries. Our method leads to a rather small number of numerically stable 
integrations. 

We apply our technique to all singular integral topologies which appear in the evalua- 
tion of NNLO double real radiation corrections to production processes of a massive system 
at hadron colliders. We present suitable phase-space parameterizations, analyze the sin- 
gularity structure of the matrix-elements, and demonstrate how to obtain their expansion 
in e with simple changes of integration variables. We then demonstrate how our technique 
can be applied to very complicated and maximally singular two- loop master integrals. In 
massless two-loop boxes overlapping singularities are very hard to treat with non-linear 
transformations, and we have not been able to find suitable ones which factorize them 
completely. On the other hand, a hybrid approach of non-linear transformations combined 
with sector decomposition is straightforward and more efficient than applying only sector 
decomposition. 

In Section 2, we review existing methods for the Laurent expansion in the dimensional 
regulator of integrals in higher perturbative orders. In Section 3 we present our method 
and we demonstrate it on typical examples of integrals with overlapping singularities in 
Section 4. In Section 5 we discuss phase-space parameterizations and the singularity struc- 
ture of double real radiation at NNLO. In Section 6 we present the numerical evaluation of 
integrals from all topologies which appear in double real radiation corrections at NNLO. 
In Section 7 we apply our method to maximally singular two-loop integrals, the crossed- 
triangle and the crossed-box. Finally, we present our conclusions in Section 8. 

2. Laurent expansion of Feynman integrals in the dimensional regulator 

Loop integrals and phase-space integrals for the calculation of production rates of phys- 
ical processes in quantum field theory are divergent in four space-time dimensions. In 
dimensional regularization, d = 4 — 2e, all divergent integrals are computed as a Laurent 
expansion in the dimensional regulator e. This task is tedious due to physical singularities, 
corresponding to infrared and collinear configurations of real and virtual particles. Sin- 
gular manifolds in the integration volume are of increasing complexity at higher orders in 
perturbation theory. We shall consider examples of physical loop and phase-space integrals 
in later sections of this paper. Here, we shall present illustrative mathematical examples 
with similar singular behavior as in realistic cases. 
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The easiest category of singular integrals is when divergences in the integrand occur 
as poles of a single variable. Consider an integral 



J o 



N f{x\,...,X N ) 



(2.1) 



with f(x) being a finite function in the integration volume. This integral is divergent for 
e = due to the pole in x\ = 0. To expand in e we use a subtraction technique, isolating 
the pole contribution, 
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In the first term, we are allowed to perform a Taylor expansion in e, given that the integrand 
is finite in the limit x± — > 0. In the second term, we can perform the integration in e easily. 
We then have 
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d N x 



-f(xi,X2, xn) ~ /(0, x 2 , ■ ■ • , x N ) 



Xl 



,n=0 



-f - / d xf(0,x 2 , ■ ■ -,x N ). 
e Jo 
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Equivalently, we write 



d N xx l 1+t f(xi, . . .,x N ) 



n! 



riog n (xi)- 









f(x 1 , . . .,X N ). 



(2.4) 



All integral coefficients in the Laurent series of the last expression can be evaluated nu- 
merically. In case of many factorized singularities, 



= I' d N x- 
Jo x- 



1 



(2.5) 



we can apply readily the same procedure, and obtain a Laurent expansion in e with the 
substitution 



die ' n\ 

n=0 



log n (xi) 



Xj. 



(2.6) 



J + 



We note that one may also encounter singularities due to poles of second or higher order, 
as for example in 
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The subtraction method can be also applied here, writing 
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(2.7) 



(2.8) 
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The integrals in x of the above expression can be computed numerically (as an expansion 
in e). 

The extraction of divergences is more complicated for integrals with overlapping sin- 
gularities. Consider as an example the integral 

/= / d Xl dx 2 - -2—. (2.9) 

JO (Xl + CX2) 

For e = the integrand becomes divergent when both xi,X2 — > 0. Here, it is easy to 
perform successively both integrations, finding the explicit result 

i r_ e - 1 -« + (i+c)--i\ (2 . 10) 



e(l + e) 

2.1 The differential equation method 

Analogous problems in realistic NNLO calculations are very hard to tackle with direct 
analytic integrations. A powerful method which has found numerous applications is the 
method of differential equations. In this approach we find a physical parameter for the 
integral and formulate a differential equation using integration by parts. In our example, 
we can write a differential equation with respect to the parameter c, by integrating the 
total derivative 

f d Xl dx 2 ^- X2 (2.11) 

J 0X2 (X\ + CX2) Z+t 



This yields the differential equation, 



91 

^~Qq ^ — -^simpler (^-12) 



Simpler = J dx 1 j- | ^ A2+e . (2.13) 



The inhomogeneous term on the right side of the above equation is simpler than /. Specif- 
ically, 

1 

{c+~x 1 ) 2+e ' 
and we find 

/simpler = ^ [(1 + C)- 1 ^ - C- 1 ^] . (2.14) 

The general solution of Eq. 2.12 involves integrals over one variable only, 



Const. + J dc /simpler (c) 



(2.15) 



thus bypassing the problem of overlapping singularities. The constant of integration can 
be determined from knowing the integral at a special value of c or by exploiting a known 
limiting behavior or scaling. For example, in our case, we could be using that 

/(l/ C ) = c 2+e /(c), (2.16) 

which we can easily derive with a change of variables x\ -H- X2 in Eq. 2.9. 
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2.2 The Mellin-Barnes representation method 

Mellin-Barnes representations allow a straightforward Laurent expansion of Feynman inte- 
grals by using Cauchy's theorem. Such representations are obtained by using the identity, 

F(N) 1 r w o+ io ° 



(A + B) N 2m 



— / dwT(-w)T(N + w)A w B~ N - w , (2.17) 



where the contour of integration is a straight line parallel to the imaginary axis, crossing 
the real axis at a point wq such that the real part of the arguments of the gamma functions 
are positive. Using Eq. 2.17 and integrating x\ and x 2 for the toy example of Eq. 2.9, we 
obtain the Mellin-Barnes representation 

i = — f*^tsm + ! ± ") Eg+ ^Itl - * -y (2.i8) 

2m J T(2 + e) T(2 + w)T(-e - w) ' v 1 

where the representation is valid (all Gamma functions have arguments with positive real 
part) if we choose, for example, e = —0.9 and Rett; = —0.2. Notice that we cannot find any 
value of Hew which renders the integral well defined if we choose e = 0. This means that 
the integral develops a pole in e. A Laurent expansion can be achieved with an analytic 
continuation method, moving the value of e from a value that the integral is well defined, 
e = —0.9 in our example, to e = and isolating with Cauchy's theorem the poles which arise 
when the arguments of e dependent Gamma functions become zero or negative integers. In 
our example, we find that T(— 1 — e — w) = T(0) develops a pole as e = —0.8 (and w = wq). 
No other pole is encountered by continuing the value of e further to e = 0. We can then 
write, 

/ = TaylorExpand(/) e=0 + Res w= _i_ e (/). (2.19) 
2.3 The subtraction method 

The differential equation method and the Mellin-Barnes method bypass the problem of 
overlapping singularities by integrating out Feynman parameters and phase-space variables 
and generating equivalent representations where overlapping singularities cannot occur. 
Both methods rely on the integration volume being well known and free of parameters, other 
than the space-time dimension. This is the case for loop integrals and inclusive phase-space 
integrations. An important class of phase-space integrals requires parametric boundaries 
which are determined according to varied selection criteria for the experimentally measured 
observables. For such integrals the differential equation and Mellin-Barnes methods are not 
generally suitable. 

One approach is to use a subtraction method in order to map the problem of fully 
differential phase-space integrations onto a problem of fully inclusive phase-space. Consider 
the toy example, 

J{XI,X2 ) 

where the function J(x\,X2) plays the role of selecting an arbitrary subregion of the in- 
tegration volume according to, for example, the wishes of the experimentalists. Using 



I[J]=j\ Xl dx 2 ( ^^ +£ , (2.20) 
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subtraction, we can re-write 

I[J} = f 1 d Xl dx 2 J( 7' Z j ~ ffi 0) + J(0, 0) [ 1 dx 1 dx 2 . 1 (2.21) 

The first integral contains only an integrable singularity as e — > and can be computed 
numerically. The second integral is a "fully inclusive" integral and may be computable 
with the Mellin-Barnes or differential equation method. 

2.4 The sector decomposition method 

A different approach is to use sector decomposition. We divide the integration region 
according to the relative magnitude of the integration variables which are required for 
the singular limit (in our example x\ = x 2 = 0), by multiplying the integrand with an 
appropriate unity, 

l = 9(x 2 >a; 1 ) + e(a; 1 >a;2). (2.22) 

This gives rise to two integration domains (sectors). In the sector with x\ > x 2 we rescale 
x 2 = x 2 xi, and in the sector with x 2 > x\ we transform x\ = x\x 2 . We then obtain 

. 1 _ e J(xi,six 2 ) , f 1 , , - 1 - e J(x 1 X 2 ,X 2 ) 



I[J} = J o d Xl dx 2Xl e (1 + ' X2)2+£ + J o d Xl dx 2 x 2 £ (1 + xi ' )2+e - (2.23) 

The singularities in both integrands are now factorized and a Laurent expansion can be 
easily achieved with a simple subtraction. The method of sector decomposition is suited 
for all types of loop and phase-space integrals. 

It is instructive to see how the method is used on a physical example. Let us consider 
the one-loop box scalar integral, 



d d k 1 



171-5 k 2 (k + Pl ) 2 (k + pi + P2 ) 2 (k + Pl + P2 + P3 ) 2 ' 

The corresponding Feynman parameterization reads, 



(2.24) 



Jo 



dx\ . . . dx^5 (1 — x\ — ... — X4) f(xi, . . . , X4) (2.25) 

Jo 

with 



/(xi, . . . , x 4 ) = \ 7 -2+ e . (2.26) 



T(2 + e) 

-SX1X3 — tx 2 Xi — ioy 

To avoid creating poles at the upper limit of the Xi integrations we apply first the method 
of primary-sectors [21]. We write 



/ dx f (]ldx i )s(l-^x i )f({x i })^26(x i -x)Y[Q(x i >x j ) (2.27) 



We now rescale 

x k = y k x, (2.28) 
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and perform the x integration. This yields 

1 / \ 2e 
I = T(2 + e) f dlft...dy 4 (X>] -. ~ ^ , nl2+e • (2.29) 

All terms in the sum can be computed in exactly the same fashion. For convenience, 
although not necessary, we use the special symmetry of this problem, y\ o y 3 and y 2 -H- 3/4, 
and cast the integral as 

/ = 2 r(2 + e) f d Vl dy 2 dy 3 (1 + yi + y 2 + y 3 ) 2e 

JO 

x { [- ayi - ty 2 y 3 ]~ 2 - e + [— tyi - sy 2 y3]~ 2 ~ £ } (2.30) 

We observe that the integral becomes singular in the following instances 

yi = and (y 2 = or m = 0). (2.31) 

We now apply sector decomposition to factorize the entangled singularity structure. We 
multiply the integrand with 

1 = o(y 2 < yi ) + e(yi < y 2 ) (&(yi < y 2 y 3 ) + 0(y 2 2/3 < 2/1)) 
= 0(2/2 < 2/1) + 0(2/1 < 2/22/3) + 0(2/22/3 < yi < y 2 ) (2.32) 

In each of the three sectors of the above equation we rescale the smallest variables with 
respect to the large ones, mapping the boundaries of the sectors to the unit cube. Specifi- 
cally, 



0(2/1 < 2/22/3) 
0(2/22/3 < 2/1 < 2/2) 
0(2/2 < 2/i) 



2/1 -> 2/12/22/3 (2.33) 

2/1 -> 2/12/2 and y 3 ->■ y 3 yi (2.34) 

2/2 -> 2/22/1 (2-35) 
We then obtain a representation of the one-loop box, 

I = h + h + h, (2.36) 
with a simple, factorized, singularity structure: 

h = 2 r(2 + e) / d yi dy 2 dy 3 (1 + yiy 2 y 3 + 2/2 + 2/3) 2e 
Jo 

x {[- syi - tr 2 ^ + [-t yi - s]- 2 -*} (mm)- 1 -' , (2.37) 

I 2 = 2 T(2 + e) / d yi dy 2 dy 3 (1 + yiy 2 + 2/2 + 2/32/i) 2e 
Jo 

x { [-5 - ty 3 ]- 2 " e + [-t - s2/3]~ 2 ~ £ } (2/22/i)~ 1_e , (2-38) 

h = 2 r(2 + e) / dyidy 2 dm (1 + 2/1 + 2/22/1 + 2/3) 2e 
Jo 

■a " i2/22/3]" 2 " e + [-* - sy 2 y 3 ]- 2 ' e } y^~ e . (2.39) 



The resulting integrals I±,I 2 , 1 3 of sector decomposition can all be expanded in e with the 
subtraction method. 
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3. Factorization of singularities with non-linear transformations 



In this section, we propose a new method for the factorization of overlapping singularities. 
We consider again the very simple case of an overlapping singularity in the two-dimensional 
toy example integral 

1 

' (c 1 x 1 + c 2 x 2 y 

We shall perform a rescaling transformation over the entire integration region, 



Ii = J dxidx 2 - ^ ^+7 



x 2 = Axi, (3.2) 

This yields, 

h = f dx lX ^ l ~ e [ X1 d\- ~~ rvoj - • (3.3) 

Jo Jo (ci + c 2 A) 2 + e 

We notice that there is a factorized singularity at x\ = 0. In this singular point of the 

x\ integration, the variable A ranges up to +oo. However, the A integration is convergent 

at e = 0, since the integrand scales as ^ for very large values of A. Therefore, we could 

immediately treat the singularity at X\ = with the subtraction method 



f'l f'OO -y 

h = dxix7 1_e / dX- — rr. . 

Jo Jo (ci + c 2 A) 2 + e [ X! 



6(A< — )-l 



1 f°° 1 

" / dX l ■ ( 3 - 4 ) 

ej (ci + c 2 A) 2 + e v ' 



We can then evaluate numerically the integrals which are produced after we perform an 
expansion in e. 

Alternatively, we could perform a transformation 1 on A to bring the integration region 
back to [0, 1]. 

X = g(u), (3.5) 
Such a transformation maps the integration region to, 

\ dx\dx 2 = dxx dug'{u) (3-6) 

Jo Jo Jg^io) 

It is very important to select carefully this transformation. A linear mapping 

g(u) = ^, (3.7) 
u 

is clearly ineffective, since it undoes the original rescaling of x 2 = Xx\. However, non-linear 
mappings, such as 

u 

9{U > X) = x + 5(l-u) (3 - 8) 
u 

9{ ^ x) = x + s 1 (i-uy 2 (3 - 9) 

x) = /ru/r fT =s ~ 1 (3 - 10) 
^(1 + 2x)(l - u) + x l 



1 Integrating A numerically, with Monte Carlo methods, requires a transformation as well, in practice, 
since one needs to generate A from some random variable that is produced in [0, 1]. 
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are effective. For almost all practical applications in this paper we employ the mapping 



u 

9 ^ x) = x + 5(i- u y (3 - n) 



with 8 often chosen equal to 1. 
Explicitly, the transformation 



xi4 _ {I + Xl ){l - x' 2 ) 

X2 ~ Xl + (i- x ' 2 y 1 ~ X2 ~ Xl + (i-x' 2 ) ' (3 - 12) 



with a Jacobian 

dx 2 Xi(l + Xl) 



dx' 2 + 

disentangles the overlapping singularity, transforming the integral of eq. 3.1 as 



(3.13) 



h= f dxidx^ 1 e (l + xi)(l-X2 + xi) e [ci(l-X2 + xi) + c 2 22r 2 ~ e (3.14) 
Jo 

The singularity in the limit x\ = and e = can be subtracted away, and a Laurent series 
expansion around e = is achieved using the expansion 



e ^ n\ 

n=0 



ln n (x) 



x 



(3.15) 



In this approach, we have achieved to factorize the overlapping singularity with a 
simple transformation. In comparison, a factorization with sector decomposition doubles 
the number of integrals, as we have seen in the previous section. Economizing in the 
number of integrals is even more significant for physical applications where entanglement 
of singularities with more variables may take place. 

Let us now revisit the one-loop box calculation using the new method instead of sec- 
tor decomposition. Following the "analytical-transformation" approach, we perform the 
change of variables on the integral of Eq. 2.30 

Vi ->• i : • (3.16) 

This yields the integral 

I = 2T(2 + e) / dy 1 dy 2 dy 3 (y 2 y3)~ 1 ~ e (l - Vi + 2/22/3)"' [2/12/22/3 + (1 + 2/2 + 2/3)(l ~ Vi + 2/22/3)] 2e 
J o 

x { [- m - t(l - j/i) - ty 2 y 3 ]- 2 - e + \-t yi - s(l - yi) - sy 2 yz]~ 2 ~ e } (3.17) 

In this integral, the singularities have nicely factorized in the term (2/22/3) -1 ~ e - I n compar- 
ison to sector decomposition, we now have to perform one integration rather than three. 
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4. Characteristic forms of entangled singularities 

In this section, we present some typical examples of integrals with overlapping singularities 
and the mappings that we use to disentangle them. 
Our second example is the integral 

h = l dxdydZ (x + U^ 
which is a trivial extension of eq. 3.1. We use the mapping 

(x + yz) : x^-^- (4.2) 

1 — x + yz 

where we have also designated the singularity structure of the integral. This mapping leads 
to 

= f 1 dxdydz (1 - x + yzf 
2 J (yz)i+* (l + yz)^ [ ^ 

where the singularities are factorized in terms of y and z. 
Next, let's consider 

f 1 1 

h= dxdydz (4.4) 
Jo (x + y + zy i+e 

Here we use the simultaneous double mapping 



which leads to 



{x + y + z) : y^ VX , z -> ZX (4.5) 

1 — y + x 1 — z + x 



f 1 dx dy dz (1 + x) 2 (l - y + x) l+e {l -z + x) 1+e 
3 " J x l +* ' ' {{l + xf-zyY+t ' { " ' 



Next, let's consider the integral 



h= dxdydzdw- — — — — (4.7) 
Jo {x + y{z + w)Y +t 



Here we use successively 



zw xyw , , . 

(x + y{z + w)) : — , x -> y — . 4.8 

1 — z + w 1 — x + yw 



The integral then becomes 

i*=f o 1 ^ 1 ^ ;v \/ 1 — v ; l2 +; 1 ~ 7 ■ (4.9) 



rl dx dy dz dw (1 + w)(l + yw)(l — x + yw) € (l — z + w ) 
y 1+e [1 - zx + w(l + y + yw) }2+t 



It is maybe instructive to see how this integral of eq. 4.7 behaves under simple rescaling. 
Consider z = X z w and then x = \ x yw. We get 

f 1 dv dw f 1 /™ [ x l yw 1 
= / <^^_ / dA / dX J: (4. 10) 
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Only the integral over X x extends to infinity and in that limit the behavior of the integrand 
is dX x /X x which vanishes at infinity. 
Next let's consider 

h= f dxdydzdw , ( xyzvj y (4.11) 
Jo (x + y + zw) 6 

Here we use the successive mappings 

/ \ xzw yzw , , „„. 

{x + y + zw) : x->- ■ , y -> ■ (4.12) 

1 — x + zw 1 — y + zw 

which brings the integral to the factorized form 

'^^. .y; (413) 







When overlapping singularities appear together with factorized singularities in the 
same variable, a slight complication appears. Consider the integral 

h= C dxdy (4.14) 



o x(x + y) 

It has a factorized singularity at x = and an overlapping singularity at x = = y. Let us 
call the singularity at x = active and the one at y = passive. In order to disentangle the 
singularity we would like to use the same mapping as in the previous examples. It turns 
out that we can do this, but only as long as we choose to remap the active singularity, i.e. 

x(x + y) : x ->• — . (4.15) 

1 — x + y 

The integral then becomes 

h= f dxdyx- 1+e y- 1+2e {l-x + y y e . (4.16) 
Jo 

which can be subtracted easily. 

Note that applying the wrong rescaling y = Xx one gets 

h= L"L (4 - l7) 

We can see immediately that the dX integral is logarithmically divergent at the active 
singularity x — > due to the upper limit of the integration region. On the contrary, 
applying the correct rescaling x = Xy one gets 



fdy! 
Jo Jo 



1/2/ (xy 

■ V»A(1 + A) ( " 8) 



which, as y — > 0, behaves as dX/X 2 which is finite. 

We see that the simple A-rescaling works as a guideline, showing when a mapping prop- 
erly factorizes the singularities of an integral. It will be instrumental in more complicated 
cases presented below. 
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Let us now consider the integral 2 

I 7 = C dxdydz ^ VZ) " - . (4.19) 
Jo xy(xy + z) 

Identifying x = = z and y = = z as two independent overlapping singularities, where z 
is passive and both x and y are active, we know from the previous example that we should 
not map z from the previous example. So one is left to map x — > z or y — > z. By the 
symmetry of the integrand it does not matter which one of these one can choose. Let us 
choose x and a slightly modified mapping that keeps the expressions simpler: 

XZ 

xy{y + z) : x -> - ■ . (4.20) 

1 — X + zx 

The integrand then becomes 

f 1 , , , (xy 2 zY(l — x + zx)~ e 

W.^ J»U + *) + (!-*)) • (4 ' 21) 

We have now managed to "activate" the singularity at z = 0. At the same time the 
singularities at y = and at x = have remained active as well. However, we now find 
an overlapping singularity at z = = y and x = 1, where the singularity at x = 1 is 
passive. Notice that we started with two independent (partially interfering) overlapping 
singularities, have treated one of them and are now left with only one, which lies at a 
different point. We shall remap z and y as follows 

y(l — x) z(l — x) 

V ^ l-y + (l-x)y ' Z ^ l-z + (\-x)z ( ' 

The integrand becomes 

Note that the remaining singularity of the integrand is integrable. 

Let's now explore the potential of a slightly different kind of mapping. We have the 
integral 

II<*VTiF (4 - 24) 

with a, b independent of x but potentially dependent on y; L . In the latter case the integral 
might have overlapping singularities, as a, b — > or x, b — > 0. We employ 



xb/a 



1 — x + b/a 

and get 

r- Prr*,* 

b N - 1 (a + b) N - 1 



(4.25) 



n ds(a(l-x) + ft)"-* 



2 We find similar singularity structures in double real radiation. 
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If N > 2 this mapping factorizes the singularity at b — > and, at the same time, exposes 
the a + b structure of the overlapping a, b — > singularity, making it ready for further 
mappings. 

Let's see, as an example, 



f 1 dx 1 dx 2 dx 3 dx< l dx 5 

+ 8 = 7 no. • \4.z<) 

Jq [X\ + x 2 x 3 + x 2 x 4 + X4X5 ) d+e 



We map: 



to get 



Xl(x 2 X 3 + X 2 X 4 + X4X5) 

xi ->■ ^ r- (4.28) 

1 - Xl + (X 2 X 3 + X 2 X 4 + X4X5) 



T I A A A A A (1 - Xl + X 2 X 3 + X 2 X 4 + X 4 X 5 ) 1+e 

h = / dxi dx 2 dx 3 dx 4 dx 5 7 — ■ — ■ ■ r^— . (4.29) 

JO (x 2 (x 3 + X 4 ) + X4X 5 ) z + e (l + X 2 X 3 + X 2 X 4 + X4X 5 ) /+e 

We can now use the mapping of eq. 4.25 with a = X3 + X4 and b = X4X5 to get 

I 8 = / dx 1 dx 2 dx 3 dx i dx 5 t— t— — — (4.30) 

Jo X4 x 5 ^ (x 3 + x 4 + X 4 X 5 ) i+e 



where F(xi) is a finite function of Xj. Noting that X4 is an active singularity, we use the 
mapping of eq. 4.25 again with a = 1 + X5 and b = X3 to get 

h = dxidx 2 dxsdx^dx 5 1+£ 1+ * 1+e . (4.31) 
J x^ x^ x^ 

Let us now see some examples where we employ a hybrid method of one-step of sector 
decomposition and nondinear transformations to factorize overlapping singularities. A 
similar singularity structure appears in twodoop massless box integrals. 

We consider 



f 1 a. 1 j ti.i 2 <'•' 3 <'•< I o 

9 = y 1 , zjz. , ~ , „.„.„_m3+ £ - ( l - ,2) 



dxi dx 2 dx3 dx4 dx$ 

[X1X3 + XiX 2 + X 2 (X 4 + X 5 + X3X4X5)] 

We split this integral in two sectors 

f 1 1 dx\ dx 2 dx3 dx4 dx5 ,.„.> 
x 2 > x 3 : / -^j- (4.33) 

JO X 2 [xi(l + X 3 ) + X4 + X5 + X3X4X5] 

which has the singularity structure of eq. 4.4 and can be factorized by the mapping of 
eq. 4.5, and 



1 1 dxi dx 2 dx3 dx4 dx$ 

x\ +t [xi(l + X 2 ) + X 2 (X4 + X5 + X4X<i)] C 



1 j. u.j-1 u,a-2 u-o-3 u,x & 

X3>X2 : ' >[ /-1 , — ^ — ; — , , m3+^ ( 4 - 34 ) 



which is of the type of eq. 4.7 and can be factorized with the mapping of eq. 4.8. 
We now move to the most complicated example of this section, the integral 



dxidx^dxi dr\ dr 2 

[x\B + X1X3A + X4T1C + X/±T 2 D + X3riT 2 S] 3 + 2e 
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with A, E finite and B,C,D finite functions of ri 2- 

We split this integral in two sectors, xi,x 4 , and we get 



xi > x 4 : I W A = J YZTTB~, 7Z ^X^T^^^'J C\ tp^ 13+^ ( L:5G ) 



dx\dx^dx^dTidT2 x\ 
[ Xl (B + Ax 3 + Cx 4 r i + L>x 4 r 2 ) + £x 3 TiT2] 3 + 2e 



which is of the type of eq. 4.1 and can be immediately factorized with the mapping eq. 4.2. 
The other sector is 

x > x . j — f 1 dx\dx- i dxA t dT\dT2 x 4 

Jq [BX1X4 + AX1X4X3 + Cx 4 Ti + DX4T 2 + £^X3riT 2 ] 3+2,: ' 

This should be further split in X3 , x 4 to get 

ldxsdX4,dTidT2 1 



x 4 > X3 : Iiobi 



f 1 dxn 
Jo 



x\ +2e [Bxi + Ax lX4 x 3 + Cn + Dt 2 + £x 3 TiT 2 ] 3 + 2e 

(4.38) 

which is of the type of eq. 4.4 and we can use the mapping eq. 4.5 to factorize it. The other 
sector is 



x 3 > x 4 : I 10B 2 



1 dx\dxzdx^dT\dT2 x 4 



x 



l+2e [Bx\XA + AX1X4X3 + Cx 4 Tl + DX4T2 + ET\T2] 3+2e 



(4.39) 

and requires further splitting. We choose to split in the variables x 4 ,ti to get 

idx3dxidridT2 1 

x 4 > ri : hoB2A 



f 1 dx\da 
Jo x\ 



, ? 2e x i+ 2e [B' Xl + A Xl x 3 + Cx 4 n + D't 2 + £riT 2 ] 3 + 2e 

(4.40) 



which is of the type of eq. 4.11, and 

" l dx\dxsdx4dTidT2 x 4 



ri > x 4 : hoB2B 



lo x\ +2e Tl +2t [Bx x x A + ^xix 4 x 3 + Cx 4 n + Dx 4 t 2 + Et 2 ] 3+2€ 

(4.41) 

which requires a final split in x 4 , T2, 

, .... ...... .... . • . \ (lT2 X 4 

T2 > X 4 : 110B2B1 



f 1 dxidx^dx^dric 



\ +2t [BxiXi + Axix 4 x 3 + Cx 4 ri + L>x 4 r 2 + Ef +2e 

(4.42) 
which is finite and 

, ^,idX4dTldT2 

x 4 > r 2 : hoB2B2 



f 1 dx\dx^d. 

Jo ~4^[- 



i +2e x\ +2e [B Xl + AE1X3 + Cn + Dx 4 T 2 + Et 2 } 3 + 2 ' 

(4.43) 

which is of the type of eq. 4.7. The original integral can be written in terms of its five 
sectors as 

ho = hoA + hoBi + IlOB2A + I\0B2B1 + hoB2B2- (4.44) 
Finally, let's consider 

I' 1 dx\dx2dxi±dT\dT2 x\ +e , . 

11 Jq [X\A + X 1 X 2 B 1 B 2 + X2X 4 Tl5 2 C + X2X4T2B1D + X2TlT 2 -E'] 3+2e 
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with A,E,C,D finite and B\2 = 2 — n^, also finite for all values of t\^. We begin by 
splitting the integral in t±,t 2 . We get 

T\ > t 2 : hiA — " ; 



J 



o [xiA + X1X2-B1S12 + X2X4T1B12C + X2X4T2T1B1D + £2Tf r 2 £ , ] 3+2e 

(4.46) 

and 



Jo 



, dx\dx2dx/±dT\dT2 x 2 

T2 > T\ : Ins - ' 



1+e 



[xi-A + xix 2 B 2 Bi2 + x 2 X4TiT2.BiC + x 2 x 4 ti-Bi 2 -D + x 2 rir|i?] 3+2e 

(4.47) 

where -B12 = 2 — t\t 2 - We notice that we can get I\\b from /iia if we exchange C and D 
and rename the dummy integration variables t\ o t 2 , so that 

WC,L>) = WAC). (4.48) 

We will use the same decompositions and mappings to factorize I\\b and I\\a (with t\ and 
T2 interchanged), so we only describe the latter below. We split I\\a in two sectors with 
respect to t\ and X4. We get 

f 1 dx\dx2dx/±dT\dT2 x 2 +e Ti 

~\ > •'' 1 : hlAi - ! ~ ~ ~ 



,0 [xi(A + x 2 B 1 B 12 ) + x 2 t?(x 4 B 12 C + x 4 t 2 B 1 D + r 2 £)] 3+2e ' 

(4.49) 

We perform the mappings 

r^Y^— (4-50) 

1 - T 2 + X4 

and then 

xir 1 2 r 2 x 4 , , 

xi -)• 5 (4.51) 

1 — xi + r 1 z r 2 X4 

which factorize all singularities. We also have 



x 4 > n : /11A2 = / 
•/ 



dx\dx2dx^dT\dT2 

2^D + nra^)! 3 ^ 

(4.52) 



[xi(^ + x 2 B[B' 12 ) + xsxIn^C + t 2j B;D + tit 2 £)] 3 + 2 ' 
We will split with respect to ti,t 2 to get 

Ti > T 2 : I\\A2a ~ 



1 dx\dX2dx^dT\dT2 X2 +£ X4Tl 

[ Xl (A + x 2 B[B'{ 2 ) + X2x\t 1 {B'1 2 C + T 2 riBiD + T 2 T 2 £)] 3+2e 

(4.53) 

which can be factorized by 

Xl -> XlX2 "' ri 2 (4.54) 
1 — xi + x 2 x|ri 

and 

Z" 1 dx\dX2dX/±dT\dT2 X 2 +6 X4T2 

T2>n : 1^26 = ^ - X2j g/// B ///) + X2X 2 TiT2(S /// C + T2B ,,, D + ^P)]3+2l 

(4.55) 
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(4.56) 



The original integral can, therefore, be factorized in six different integrals: 



hi = hlA + hlB = hlAl + hlA2a + hlA2b + {u t) 



(4.57) 



5. Double real radiation for final states with massive particles 

One of the major challenges at NNLO in QCD has been the computation of the double real 
emission part of the cross-section. While the computation of the matrix elements with iV+2 
particles in the final state is not a problem per se, difficulties arise when one integrates over 
the phase space of the two unresolved particles. The corresponding integrals are infrared 
divergent in the soft and collinear limits and are dimensionally regulated. The divergences 
have to be subtracted before the integrals can be numerically evaluated. As long as the 
singularities are factorized, as they usually are at NLO, it is straightforward to use a 
Laurent expansion over the singular variables, and evaluate its coefficients numerically. At 
NNLO, the singularity structure of the integral is more intricate, as line and overlapping 
singularities appear, and the desired factorization is not straightforward. 

The method of sector decomposition has already been applied successfully to achieve 
this factorization for hadron collider [7, 8] and decay processes [24-26] . A drawback of the 
method is that it leads to a large number of sectors. The goal of this paper is to replace 
sector decomposition for double-real radiation integrals with an economical factorization 
method based on non-linear transformations. 

5.1 Infrared singularities in double real radiation 

We consider double real emission to a generic NNLO 2 — > n + 2 process (see Fig 1) with 
n massive particles and 2 massless partons in the final state. We denote the momenta of 
the incoming particles by qi,q2, those of the outgoing massive particles by p\ . .p n and those 
of the two unresolved partons by and q^. Infrared singularities in this phase space will 



occur whenever q% and/or q^ become soft or collinear to qi, q<i or to each other. For the 




Pn 



Figure 1: Notational setup: qi are massless particles, while pi are massive. 
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case of double real radiation to the production of a single massive particle (e.g. Higgs, W 
or Z production) potentially singular propagators can be summarized as 



S34 = 2q 3 .q± 

513 = -2q 1 .q 3 

523 = -2<7 2 -<?3 

514 = -2q 1 .q 4 

524 = -2q 2 -q4 (5.1) 

and 

S134 = + 94 - qi) 2 = S34 + S13 + S14 

S234 = q2) 2 = S34 + S23 + S 2 A- (5.2) 

Note that S123, S124 are bounded from below. Further soft singularities can be found if 
there are colored massive particles in the final state, which can radiate off soft gluons. One 
can then get also the following possibly singular denominators: 

hi = 2q 3 .pi 

tu = 2q 4 .pi (5.3) 

and 

*34i = (<?3 + 94 + Pk) 2 ~m\ = S34 + t 3i + Ui (5.4) 

for i > 1. Since 

hi = 2q 3 -Pi = 2E 3 (Ei - |pj cos6» 3i ), > |pj (5.5) 

the soft singularity is factorized in E 3 . Whenever some heavy colored state radiates off 
two gluons we can also get the denominator t 3 a, it can only become singular in the double 
soft limit when E 3 = = £4. However the double soft limit will always be factorized as 
we will show in the next section. 

Let us now discuss the denominator structure of the most singular diagrams which 
one could expect in double real radiation: those where radiation is emitted by initial state 
particles. We have illustrated the propagator structure of these topologies using some 
diagrams containing gluons in Fig. 2 (diagrams containing just massless quarks correspond 
to the same topologies). 

Diagrams whose propagator structure can be related to the ones in Fig. 2 by a simple 
interchange of q 3 with q^ or of q\ with q 2 will also fall into the same topology. 

By considering square and interference terms of the topologies C\ , C2 and C 3 , we obtain 
the following list of integrals: 



1. Topology Ci (8) C\\ 



r d<s> 3 N({ StJ }) ^ r d<f> 3 N({ Sij }) 
J (S13S24) 2 J S13S23S14S24 



(5.6) 



Ci : 




Pl.,7. 



1/(S13S24) 



92 wmrmm' 94 



'h 



Pl..n 



C 3 : 




l/(s M s 



234 J 



1/(«134«13) 



l/(*31«4l) 



Figure 2: Most singular topologies:Ci, Ca, C3, C4 



2. Topology C 2 <g> C 2 : 



3. Topology C 3 ®C 3 : 



(S34S134) 2 ' J s| 4 Si34S 2 34 



d$ 3 N({ Sij }) f d$ 3 N({ Sij }) f d$ 3 N({ Sij }) 



r d$ 3 N({ 8ij }) ^ r 
J (S13S134) 2 j 



Sl3-S23'Sl34'S234 J S13S24S134S234 



4. Topology C\ <g> C 2 



5. Topology C\ (g) C 3 : 



6. Topology C 2 <8> C 3 



I 



d<S> 3 N({ SlJ }) 



S34'S234Sl3'S24 

d$ 3 N({ Sij }) f d$ 3 N({ Sij }) 



Sl34'Sl3'523Sl4 J Si34S 13 Si4 

d$ 3 N({ Sij }) f d$ 3 N{{ Sij }) 



• s 34S 1 3 4 Sl3 J S34S134S234S23 



(5.7) 

(5.8) 

(5.9) 
(5.10) 
(5.11) 
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7. Topology C 4 <g> C 4 : 



8. Topology C 4 ®Ci: 



9. Topology C 4 <g> C 2 : 



10. Topology C 4 ® C 3 : 



/ 



d$ 3 #(K}) /• d$ 3 i V({«ij}) 



A 2 f 2 



/• c2$ 3 A 

7 ii3*7 



/ 



d<E> 3 iV({^}) 
ii 3 ij4Sl 3 Sl 4 

d<D 3 A({^}) 
ii 3 ij4S 3 4Sl34 

d<D 3 A({^}) 
ii 3 ij4Sl 3 Sl34 



(5.12) 



(5.13) 



(5.14) 



(5.15) 



Where d$ 3 is the differential double emission phase space element for 2 + n final state 
particles, and N({sij}) is in general a finite function of the kinematical invariants. 

The topology C 4 contains only soft singularities similar to those in C\ . The topologies 
C 4 ® C 4 and C 4 (8) Ci are, therefore, easier than C\®C\. They can be treated exactly like 
C\ C\ and we will not discuss them in what follows. 

5.2 Phase-space of double real parton radiation 

We would like to point out that different parameterizations of the phase-space can factorize 
different sets of kinematic invariants. We will derive two such parameterizations which 
allow for a more convenient numerical evaluation of diverse diagrams, according to their 
topology. 

The phase-space of n massive particles in four dimensions is: 

/ n \ n 

d$ n (>A;mi, ..,m n ) = (2vr) 4 - 3 " JJdV+Cp? - m?) 5 (4 \ gi + q 2 - (5.16) 



i=l 



where s = {q\+q2) 2 - We assume that a 2 — >■ n process exists at leading order in perturbation 
theory, and a strictly four-dimensional evaluation is therefore sufficient. At NNLO, the 
double emission phase space is given by including two further massless particles (whose 
momenta we denote by g 3 and g 4 ) 



d$ n+2 ( > /I;mi,..,m n ,0,0) = (2tt) 4 - 3 " \J[ d A Pi 5 + ( P * - m?)J ^vO^V^+feV^+tel) 

x <J (d) (gi + 92-^^-93 -94). (5.17) 



i=i 



We factorize the double real phase space into a 3-particle phase space times an n-particle 
phase space as follows 

/ds 
—^d^3(^/s; 0, 0, y/si.. n )d$n(y/ s i..n; mi, .., m„). (5.18) 
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* P2 



Figure 3: Phase space factorization 

with si.. n = (X^iLiP*) 2 shall denote the center of mass energy (or invariant mass) of the 
n massive momenta pi,..,p n . This is depicted graphically in Figure 3. The limits of 
integration of si.. n are 



mi 



(5.19) 



S > Si.. n > I X/ 

\i=l / 

and parameterizing si.. n linearly we get 

SLn= | s " |E m, j j £5+^^771^ . (5.20) 

The parameter X5 6 [0, 1] then uniquely defines the double soft limit when X5 = 1. In 
the following discussion we will use the variable 

z = 

s 



(5.21) 



which in the special case of n = 1 reduces to z = — . Then the variable X5 is trivially 
removed and the double soft singularity occurs whenever s = . 

In the following we will assume that one can parametrize the n-particle phase space 
d& n , and we will focus on the phase-space of the potentially unresolved massless partons 

5.3 Energies and angles parameterization 

The three particle phase space element d<3?3 is 

d$ 3 (V^;0,0, ^7n) = (2vr) 3 - 2 Vg 3 5 (+) fe 2 ) ( i^45 (+) (g4 2 )^Q< J (+) (Q 2 -si.. ri )<5 d (g 1 +g 2 -g3-g4-g). 

(5.22) 
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Integrating out Q and using that d d q8^ + \q 2 ) = dEE d ~ 3 dQS d -^ /2 we get 

d$ 3 (v^;0,o,VSI^) = (27r) 3 - M ^^ 1) d^ 1) ^ 3 ^4(^3^4)^ 3 

x 5 {+) {s - ai.. n - 2^(^3 + E 4 ) + 2^3^4(1 - cos #34)). (5.23) 

We can solve the delta constraint for the energies in a symmetric way using the following 
ansatz: 

E3 = 2"^( 1 ~ z ) x i K 

E A = iv^(i-z)(l-xi)K. (5.24) 



We find 



1 - VI - 2(1 - z) Xl (l - - cos^I) 

(l-z)xi(l-Xi)(l-COSe34) ) 1 J 



The double soft limit now appears when z — > 1, while the single soft singularities occur as 
a?i — >■ 0, 1. After this transformation the phase space volume becomes 

= (2vr) 3 - 2d dVLf~ 1] d^~ 1] d Xl K(l - z) f s(l - z) 2 k 2 Xi (1 - x{ t 



16^/1 - 2xi(l - xi)(l - z)(l - cos6> 34 ) V 4 

(5.26) 

Having solved the energy constraint we move on to parametrize the angles. Choosing 
the z-axis as the direction of qi, we directly parameterize the angles which q% and q% make 
with the z-axis. Finally we parametrize the angle 4> between q% and 54 in the x-y plane 
leading to the following expressions of the solid angles 

dn ( 3 d ' 1] = dnl d ~ 2) dcos9 3 {sm9 3 ) d - 4 

dn{ d ~ 1] = dn { 4 d ~ 3) dcose 4 {sm9 4 ) d - 4 dcos^{sm^) d - 5 . (5.27) 

Suppressing any extra dimensional components our 4-vectors are then fully parametrized 
as Q3 = #3(1, sin #3, 0, cos #3) and q$ = £4(1, sin 64 sin <fi, sin O4 cos (f>, cos 64). Mapping the 
remaining angles linearly, i.e. cos #3 = 2x3 — 1, cos #4 = 2x4 — 1 and (p = X2vr, one obtains 

3+2e n ^(1-z) 3 k 4 xi(1-xi) 

v ~ 2 — K 

(s 2 (l - z) Vx 2 (l - xi) 2 x 3 (l - x 3 )x 4 (l - x 4 ) sin 2 (7rx 2 ))" e (5.28) 



f (2ir) i+2e f 1 

J d ^ 3 = 16r(l - 2e) J dx ^ dx t dx ?>dx4 



The following lists the propagators of massless partons in this parameterization: 

513 = s(l - z)kxix 3 

S23 = s{l - z)kxi(1 - x 3 ) 

514 = — s(l — z)k(1 — Xl)x4 

s 24 = - 8 (1 - z )k(1 - xi)(l - x 4 ) (5.29) 
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and 



«34 = S(l - z) 2 K 2 Xi(l - Xi)x 34 
S134 = -5(1 — z)k [(1 — z)kX\(1 — Xi)x34 — X1X3 — (1 — Xi)x4] 

S234 = s(l - z)k [(1 - z)kxi(1 - xi)x 34 - xi(l - x 3 ) - (1 - xi)(l - x 4 ))] (5.30) 

where 

^34 = X3 + X4 — 2X3X4 — 2 COs(x27r) \J ' X%{1 — X3)x4(l — X4) (5.31) 



and 

1 - v / l-4(l-2;)xi(l-xi)x34 

1^ — — " 

2(1 - z)xi(l - xi)x 34 
The angle between q% and q^ is related to 



(5.32) 



1 — cos #34 1 — cos #3 cos #4 — cos <j> sin #3 sin #4 
^34 = 2 = 2 ' ^ ^ 

This expression exposes the weak point of this parameterization: it gives rise to an over- 
lapping line singularity when cf> = and #3 = 64 i.e. when qs is parallel to q^. Nevertheless 
the above construction can be used to fully subtract all phase space integrals which do not 
contain singularities in X34, i.e. which do not contain S34, 5134,3234. 

Let us now analyze the singularities in this parameterization. While si3,S23>si4 and 
S24 are fully factorized, there is a overlapping line singularity in S34 when X34 = 0. 
Furthermore there are overlapping singularities in S134 and S234. For S134 there are 3 
different possibilities 

a) X3 = and X4 = 

b) X3 = and x\ = 1 

c) X4 = and x\ = (5.34) 
while for S234 the singularities are located at 

a) X3 = 1 and X4 = 1 

b) X3 = 1 and x± = 1 

c) X4 = 1 and xi = 0. (5.35) 

We can now apply this parameterization to all integrals of type C\ <8> Ci,C 3 (8) C3 and 
Ci®C 3 . 

5.3.1 Line singularities in the energy and angles parameterization 

One can use a non-linear transformation to get rid of the overlapping structure in X34 [27]. 
A convenient way to derive such a mapping is remapping X34 from Xg 4 = xja($ = 0) to 
^34 = ^34(0 = 1) using 

~ _ X 34 X 34 /r oc\ 

^34 - Tq— ,_ + _ : (5.36) 

X 34 X2 1^X34 X 34 J 
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It is then apparent that X34 will vanish whenever Xg 4 or will, for any value of x 2 - 
And that the overlapping line singularity is then re-casted into just a line singularity. To 
aid numerical stability we perform the mapping X2 — > (1 — cos(x27r))/2, such that X34 
becomes 

X34 = , _ . (5.37) 

x 3 + x 4 - 2x3X4 + 2cos(x 2 7r)i/x 3 (l - x 3 )x 4 (l - x 4 ) 



This is in fact identical to the mapping in [28]. The phase space volume then becomes 

(2vr)- 3 + 2 - [\,,, ( s(l-z)W Xl (l- Xl ) \ 
161X1 - 2e) J dx ^ dx * dx * { ^ J 

/ \ l-2e 

(s 2 (l - z)Vx 2 (l - X!) 2 x 3 (l - x 3 )x 4 (l - x 4 ) sin 2 (7rx 2 ))" £ ( , XM , I (5.38) 

V|x 3 -x 4 |/ 



To factorize the line singularity in S34 (at X3 = X4) we are forced to split the integration 
region in two, separating X3 < X4 from X4 < X3. 

5.4 Hierarchical parameterization 

Since in the energy and angles parameterization the invariants S34, S134, S234 had line and 
overlapping singularities, it is worth having a second parameterization which factorizes 
these, but may not factorize the others. Our second parameterization closely resembles the 
features of the rapidity parameterization published in [7], however it is somewhat simpler. 
In this parameterization the three particle phase space element d^s is 

d<Mv^;0,o, VsLn) = (2tt) 3 - m ^ (+ ^ 3 V^ (+ ^ 

(5.39) 

is first factorized into a product of two 2-particle phase spaces 

d$ 3 (v / s, 0, 0, = J ^^ 2 (v / 5, Vm, V^)d$ 2 (Vm, 0, 0). (5.40) 

with 

d$ 2 (Vs, v^I, ^~ n ) = (27r) 2 - d d d Q5^\Q 2 - Sl .. n )d d Q5^\Q 2 - s u )S d {qi + q 2 - Q - Q) 

(5.41) 

and 

^2(^34,0,0) = (27r) 2 - d d%6M(q 2 )d d q 4 6^(q 2 )5 d (Q - q 3 - q A ). (5.42) 
We can parameterize cZ<&2 ^/s34, y/si.. n ) in terms of S134, yielding 

^ 2 (v^, V^34> v 7 ^) = (27r) 2 - d j-dn d ~ 2 (Q ± ) d - 4 ds lu . (5.43) 
To satisfy Q± > 0, we take 

n , , S 1U (s + S134 - SL.n) 

U < S34 < 

S134 - Sl.. n 

> si34> («i..n-a)- (5-44) 
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d^2(y/sM', 0,0) can be parameterized in terms of the invariants si 3 and S23 yielding 

1 ■ - — d— 3 r/„ A „;„ 



^2(^/534; 0,0) = (2vr) 



2-cZ 



8Q ± s 



dsi 3 ds 2 3dQ d d [(j> 3 )± sin 



(5.45) 



where is the angle between (p 3 )± and Q±. We fulfil the constraint {p 3 )± sin0 > to find 
the limits of integration for si 3 and then for S23- 
Parameterizing S134, S34, S13 and S23 linearly we arrive at 

d$3 = T^rvl iT\ / dxidx 2 dx 3 dx A -± — ' y — ' 

16r(l-2e)7 V z + x 1 (l-z) 

' ,s 2 (l - z) 4 xj(l - xi) 2 x 2 (l - X2)a; 3 (l - X3) sin 2 (7TX4 
Z + Xl(l - z) 

The invariants in this parameterization are 

s(l - Z) 2 X\{\ - X\)X2 



I 



(5.46) 



S34 
S134 
■S234 

S2A 



Z + Xl(l — z) 

-s(l — z)x\ 



-s(l — z)(l — Xl 
-s(l - z)(l - xi)x 3 



z + Xi(l - x 2 )(l - z) 



Z + Xl(l - z) 

-s(l-z)(l-xi)(l-x 3 ) 



(5.47) 



and 

si 3 = -s(l - z)zi 

Sl4 = — s(l — z)xi 



2:3(1 - £2) + 



(l-x 3 )(l-x 2 ) + 



3:2(1 - 3:3) 

Z + Xl(l - z) 

X2X3 



2 COS (7^4)1 



/ 3:2(1 ~ ^2)^3(1 ~ ^3) 
2 + Xl(l - z) 



Z + Xl(l - z) 



2cOs(7TX4)i 



1x2(1 - x 2 )x 3 (l - x 3 ) 

Z + Xl(l - z) 



We see that the only invariants which are not factorized are S13 and S14. The variable 
Si 3 contains overlapping singularities at X3 = = X2 and X3 = 1 = X2 as well as an 
overlapping line singularity at X4 = 0, X\ = 1,X3 = X2, while S14 contains overlapping 
singularities at X3 = 0, X2 = 1 and x 3 = 1, X2 = 1 as well as an overlapping line singularity 
at X4 = 1, xi = 1, X3 = 1 — X2- 

5.4.1 Line singularities in the hierarchical parameterization 

Consider the expressions 

J{Pl,Pl,PZ,P±) J(Pl,P2,P3,P4) 



S13S24 



S13S23 



(5.48) 



with J(pi,P2,P3,Pa) a finite numerator function. They both contain a line singularity due 
to sis i n the denominator. We now use the partial fractioning identities, 



1 



1 



S134 f23jt\ 

S13-524 S13S234 + S134S24 V s 13 S24, J 



(5.49) 



24 



! ^134 + £234, (5 5()) 



S13S23 Sl3«234 + S23S134 V s 13 s 23 

The term S13S234 + S134S24 has an overlapping singularity at X3 = 1 = x 2 , while the term 
S13S234 + S134S23 nas an overlapping singularity at X3 = = x 2 . Then we exchange 1 -H- 2 
and 3 o 4 in the term containing S13 to rotate the line singularity out, i.e. 

J(Pl,P2,P3,P4) = J(Pl,P2,P3,P4) + J(P2,Pl,P4,P3) S 2 34 , g ^ 

Sl3«24 S13S234 + S134S24 S24 

J{Pl,P2,P3,P4) = J(Pl,P2,P3,P4) + J{P2,Pl,P3,P4) S234 , g ^ 

Sl3^23 S13S234 + S23S134 S23 ' 

and we are left with just overlapping singularities, which can be treated as explained in 
the following section. This trick was first discovered by Frank Petriello [29] and it has been 
used in the implementation of the program FEHiP described in [7], it has been also been 
used in the evaluation of doublereal counterterms in [30]. 

6. Numerical evaluation of double-real radiation phase-space integrals 

In this section, we present a numerical evaluation of all types of scalar phase-space integrals 
which appear in NNLO calculations. To evaluate our integrals numerically we choose the 
point (s = 1, z = 0.1). We will use the notation 

Xi = l — Xi, (6.1) 

where the Xi G [0, 1] are parameters of integration. 

1. Topology Ci <g) C\\ 

(a) The integral 

/ — ^— (6.2) 

J S13S23S14S24 

fully factorizes in the energies and angles parameterization (Section 5.3), we 
obtain 

hla = 0.09400(2) + a01 °f 1(4) - - °-0°m^4946(l). (M) 

(b) The integral 

Jm= /^3 (S34S ;T 23)2 (6-4) 

J S 13 S 24 

with the numerator structure as in [7], factorizes in the energies and angles 
parameterization (Section 5.3), we get 

7, u - 0.023885(3) + °°°^0 6 (3) + 0.00036930(4) ^ 

2. Topology C 2 ®C 2 : 



hla 
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(a) The integral 

factorizes in the hierarchical parameterization (Section 5.4). The numerator 
function has the scaling behavior N({sij}) ~ S34S134 [7]. We obtain 

!«.= [ = 0-0011728(1) - °"°° 5 ° 726 « _ 0-000125982656(0) 

(b) The integral 

i=f¥W (0.8) 

factorizes in the hierarchical parameterization (Section 5.4). The numerator 
scales as N({sij}) ~ S34. We obtain 

/» = / ^ = -0.0015003(2) + 000112726 ' 5 » + O-OOO^lMO('). 

7 S34S134S234 e 

(6.9) 

3. Topology C 3 C 3 : 
(a) The integral 

J S 234 S 24 

factorizes in the hierarchical parameterization (section 5.4). The numerator 
structure can be found in [7]. We obtain 

/ 33a = -0.00384 1 (2) + a0007 e 8 ' 4 ' 4 > + a0001 t f 5 '" . (0.11) 



(b) The integral 

'336 



(6.12) 

S(L34S234'Sl3 s 23 

neither factorizes in energies and angles nor in the hierarchical parameteriza- 
tion. We use the hierarchical parameterization (section 5.4), since fewer over- 
lapping singularities are present there. Using partial fractions, as described in 
section 5.4.1 we can rewrite the integral as 



^336 



I 2* 

J S23«134(S134S23 + «234Sl3) 



This contains the following substructure 

11 

X3 x 3 A + X2X3B + 0^x2X2X3X2 

with A, B, C finite. This becomes singular when X3 = = X2 where X3 is active. 
We factorize this singularity by applying 

(1 -X3) +X 2 
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and obtain 

l» = 0.023155(3) + a0076371(1) _0-000«W61M6(1) 
(c) The integral 

fee = f ^ (6-17) 

J S134S234S13S24 

is similar to /33b in the hierarchical parameterization (Section 5.4). Partial 
fractioning as before we get 

'».= / v <« 8 > 



This contains the substructure 

1 1 



£3 x 2 ^4 + X3X2B + 0^x2X2X3X3 



(6.19) 



with A, B, C finite. This becomes singular when X3 = = X2 with X3 being 
active. We disentangle this singularity by applying 



* 3 -> n ??, - • (6-20) 
(1 - x 3 ) + x 2 



We then obtain 



I, 3C = 0.12567(9) - gl - ^ + ° »°^9%12364(0) (g ^ 



4. Topology Ci (g> C 3 : 



(a) The integral 

7 S 2 34S24 S 13 

with the numerator structure as in [7]. The singularities factorize in energies 
and angles. We immediately obtain 

l aa = -0.0040885(4) _ ° ~30 (1 ) ^ 

(b) The integral 

hs b = I d*s Ni{S - }) (6.24) 
J S134S13S14S23 

has a quadratic divergence due to the term S134S13S14. This means that N({sij}) ~ 
{■5134, S13, S14}. Such that 

/„ = {/-*•-, (6.25, 

U S 13«14S23 J S134S14S23 J 

the first of which is a sub-topology of C\ while the second is a sub-topology of 
C 3 2 - 
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5. Topology Ci <g> C 2 : 
(a) The integral 



(6 . 26) 

J S34S234S13S24 

has a quadratic divergence due to the term 5345234524- The numerator can have 
the following scalings: N({sij}) ~ {534, S234, «24}- We therefore consider the 
following possibilities 

, 12 ={/^^,/^^,/^?^.} (6.27) 

U S 34S13«24 J S34S234S13 J S234S13S24 J 

The last of these is a sub-topology of Cf and does not merit further attention. 
We will evaluate the other two in the hierarchical parameterization. For 



Il2a 



[ d * 3 (6.28) 
J S34S13S24 



we use the same strategy as we used for /33c. We obtain 

/,* = 0.07115(1) + 00069i " 5 ' 1 ' - a0 ° 29 2 912 "> - 0000839 f 709 "". (6.29) 
The integral 

/i26 = / d * 3 (6.30) 

7 S34S234S13 

factorizes in the hierarchical parameterization. We obtain 
7 126 = 0.0198554(9) + - °^8965(4) _ 0.000279961236(1) 

(6.31) 

6. Topology C 2 ®C Z : 
(a) The integral 



/*W<*» (6.32) 

J S34S234 S 23 



^23a 

S 34S234 S 23 

factorizes in the hierarchical parameterization, but carries a cubic divergence in 
(1 - xi) ~ S234- Taking N{{sij}) ~ s| 34 , we get 



I- 



23a 



[ (6.33) 

y S34S23 



which just is a sub-topology of I\2a- While other numerators are possible these 
do not give different singularity structures. 

(b) The integral 

^ = r wM) (6 . 34) 

J S34S134S234S23 

factorizes in the hierarchical parameterization, but carries a quadratic divergence 
in (1 — x\) ~ S234. A minimal choice for the numerator is N({sij}) ~ £234 in 
which case we recover 7 2 2b- Hence no new singularity structures can be obtained 
from this topology. 
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7. We will now consider interferences of C4 with Ci and C3. One can evaluate these 
interferences in the energies and angles parameterization. In the following we will 
use £13 ~ E3 ~ (S13 + S23) an d £24 ~ £4 ~ («14 + «24)- 

(a) Topology C2 £3 C4: 
The integral 

i 24 = f -, r (6.35) 

J S34S134(S13 + S23)(«14 + S24j 

has the following singularity structure 

1 1 , 

(6.36) 



XiX\(x3 — X4) Ax\Xi(x3 — X4) 2 + BX1X3 + Cx\X^ 

in the energy and angle parameterization after the mapping (see 5.3.1) is applied. 
We first split the integration region into two sectors which we define as X3 < X4 
(sector 1) and X4 < x 3 (sector 2). After this sector decomposition we are still 
left with overlapping singularities at X3 = = x\ in sector 1 and at X4 = = X\ 
in sector 2. These can be disentangled using 

- , 21X3 (R QT\ 

x\ ->■ — r (6.37) 

(l-xi) + x 3 V 1 



in sector 1 and 



in sector 2. We then obtain 



xi -»• - XlX . 4 (6.38) 

(l-xi) + x 4 V ; 



/« = -o.oo 6 % 6(3 ) - a ° 010708 ' 3 ' + "■ 0006 f°« 1 ' + o-o°o»raww4(o), 

(6.39) 

(b) Topology C 3 0C 4 : 
The integral 

% = / . f 3 , ( ~ ; (6.40) 

J S13S134(S13 + S23j(Sl4 + S24J 

has the following singularity structure 



1 1 



X1X3 ^4xiXl(x3 — X4) 2 + -BX1X3 + Cx\X4 



(6.41) 



It contains no line singularity but several overlapping ones located at X3 = = 
X4, X4 = = x\ and at x 3 = 0, X\ = 0. To separate the two singularities we first 
partial fraction the soft singularities by multiplying by 1 = xi + X\. We then 
treat the two terms with different nonlinear transformations. For the first term 
we apply the mapping 

^ £3x4x1 . . 

x 3 -> -p. r-: — ■ (6.42) 

(1 - X 3 ) + X4X1 

since X3 is the only active singularity, it is clear that we had to remap it. The 
second term is more difficult, since both X3 and x\ are now active. We apply 
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the following sequence of mappings: 
First let 



X3X4 (P. A1\ 

X3 n ( 6 - 43 

(1 - x 3 ) + x 4 



and then 



xi 



xix 3 



(1 - Xl) + x 3 



x 4 -> • (6.44) 

(l-x 4 ) + x 3 V 7 

(6.45) 

We obtain 

/„ = -0.32519(4) - a ° 48942 ' 2 ' - °° 062 f 7 ' 3 > - °^9922473(3) 

e e 2 e 3 

7. Two loop examples 

In what follows we show how one can use non-linear mappings to disentangle singularities 
in two-loop integrals appearing in NNLO virtual amplitudes. We treat only two indicative 
cases, the massless non-planar triangle with one leg off-shell and the massless non-planar 
box with all legs on-shell, due to their particularly intricate singularity structure. Integrals 
involving masses are in general simpler as far as factorization of singularities is concerned. 

7.1 The massless non-planar triangle with one leg off-shell. 

The two- loop, non-planar triangle with one off-shell leg (see fig. 4) and momenta pi,P2 

, = f d d ki d d k 2 1 

™ J m d / 2 m d l 2 kl{k l +p l fkl{k2+P2) 2 {ki+k 2 f{k 1 + k 2 +Pi+P2) 2 { '' 

A Feynman parameterization reads: 

Xtri = 4 2 + 2e I' d Xl dx 2 dzdydx r - ~ V }~ 1 ~^\~ Z) ~\l +2e (7.2) 

Jo [x(l — x) + yz(x — xi)(x — x 2 )\ 




Figure 4: The massless non-planar two-loop triangle with one legs off-shell 

The first overlapping singularity is at x = or x = 1 and y = 0. We also notice that 
there is a singularity at y = 1. To avoid infinite looping we must first guarantee, as in 
sector decomposition, that no singularities occur at the upper limit of integration. 
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We split x in the two intervals R a = [0, 1/2] and Rb = [1/2, 1] and map the integration 
region back to the unit hypercube. In R a , 

x -)• x/2, 

and in Rb, 

x — > 1 — x/2, x\ — > 1 — x±, X2 — > 1 — X2- 

This gives two identical integrals (the integral in eq.7.2 is invariant under the combined 
x — > 1 — x, X\ — > 1 — xi, X2 — )■ 1 — X2), and we can write 



Jo 



zy 1+t (l - y)- l - e {l - z)- 1 -' 
lo [ x ( 2 - x) + yz(2xi - x)(2x 2 - x) l2 " :2f 



Xtri = 4 2+2e / dx x dx 2 dzdydx ^ i ^ i '- 3+2" ( 7 - 3 ) 



Note here that the denominator of eq. 7.3 has the same singularity structure as 

x + yzx\x 2 - yzx(x\ + x 2 ) (7.4) 

In particular, there are still singularities at upper corners of the hypercube, when 
y — > 1 and z — >■ 1. This leads us to split the y integration region, y — > y/2 and y — > 1 — y/2 
and the z — > z/2 and z — > 1 — z/2. 

We obtain 

Xtri = Xtriu + Xtrn 2 + Xtri 2 i + Xtri 22 , (7.5) 

where 

Xtriu = 2 6+9e dx x dx 2 dzdydx V +£ (2 - j/)' 1 '^ - z)' 1 ^ 

Jo [4x(2-x) + yz(2x 1 -x)(2x 2 -x)] 2+2e 

Xtri 12 = 2 6+9e dx x dx 2 dzdydx (2 - % 1+e (2 - y)' ^V 1 ^ 

7o [4x(2-x) + y(2-z)(2xi-x)(2x 2 -x)] 2+2e 

Xtri 21 = 2 6+9e Z 1 dx x dx 2 dzdydx z(2 - y) 1 ^y~ l ~ f -(2 - z)^ 

7o [4x(2-x) + (2-y)z(2xi-x)(2x 2 -x)] 2+2e 

(2 - z)(2 - y^+V 1 "^ 1 ^ 



Xtri 22 = 2 6+9e / dx\dx 2 dzdydx- 
Jo I 



[4x(2 - x) + (2 - y)(2 - z)(2xi - x)(2x 2 - x)] 2+2e 

(7.9) 

The first three sectors, Xtriu, Xtri\ 2 , Xtri 2 \ are free of singularities at 1. They can 
be directly treated by a non-linear mapping. 

Xtriu has a singularity structure equivalent to that of x + yzx\x 2 , similar to eq. 4.1 
and we use the mapping (directly analogous to eq. 4.2) 

xyzxix 2 

x — > . (7-10) 

1 - x + yzx\x 2 

Xtri\ 2 has a singularity structure equivalent to that of x + yx\x 2 , also similar to eq. 4.1 
and we use the mapping 

xyx x x 2 

x ->■ . ('•!!) 

1 - x + yx\x 2 
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Xtri2\ has a singularity structure equivalent to that of x + zx\X2, also similar to eq. 4.1 
and we use the mapping 

x -> XZX1X2 . (7.12) 

1 — X + ZX1X2 

Xtr%22 is a bit more complicated. Its singularity structure is the one of 

x + xix 2 - x(xi + x 2 ) (7.13) 

It retains singularities at x± t 2 — > 1. We therefore split this integral further in xi — >■ xi/2 
and x\ — > 1 — £i/2 as well as x2 — >■ x2/2 and X2 — >■ 1 — X2/2. 
We obtain 

Xtri 22 ii = 2 5+9e f 1 d Xl dx 2 dzdydx- {2 ~ Z)(2 ' ' rr 11 ' 



jo [4x(2 - x) + (2 - y)(2 - z)( Xl - x)(x 2 - x)] 2+2t 

(7.14) 

with a singularity structure equivalent to x + X1X2 (i.e. eq. 4.1) for which we will use 
the mapping 

XX!X 2 ,_ , 

x -)• . (7-15) 

1 — X + X1X2 



(2- z)(2-y) 1+e y- 1 - e z- 1 " e 



Xfrz9 2 i2 = 2 5+9e / dx\dx2dzdvdx , ., 

Jo [4x(2-x) + (2-y)(2-z)(x!-x)(2-x 2 -x)] 2+2e 

(7.16) 

with a singularity structure equivalent to x\ + x(xi + y + z + X2), similar to eq. 4.7 for 
which we will use the following sequence of mappings 

xix yx 2 zxi xix 2 

x\ ->• 1 ■ — , y -»• ^ ■ , z -»• ■ , xi ->• ■ (7.17) 

1 — xi+x 1 — y + X2 1 — z + xi 1 — X1+X2 



Xtri222\ = 2 5+9e / dx\dx2dzdydx- 
Jo 



(2 - z)(2 - y y+*y 



l+e,.-l-e_-l-e 



[4x(2 - x) + (2 - y)(2 - z)(2 - xi - x)(x 2 - x)] 2+2e 

(7.18) 

with a singularity structure equivalent to X2 + x(x2 + y + z + xi), similar to eq. 4.7, for 
which we will use the following sequence of mappings 

x 2 x yxi zx 2 x 2 xi 

X2 -»• ^ . — , 1/ ->• 1 ; 1/ , z -»• ■ , x 2 -»• ■ (7.19 

I-X2+X l — y + xi I-Z + X2 I-X2+X1 



5+q f f 1 , , , , , (2-z)(2-y) 1+e y- 1 - e z- 1 " e 

[4x(2 - x) + (2 - y)(2 - z){2 - x x - x)(2 - x 2 - x)] 2+2 ' 



Xtri2222 = 2 5+9e / dx\dx2dzdydx 
Jo 

(7.20) 
which is finite! 

We therefore end up with 7 different integrals to be numerically evaluated. This should 
be contrasted with the 64 number of sectors one arrives using sector decomposition. The 
numerical convergence of these integrals poses no additional problems and we have checked 
that the numerical result agrees with the analytic result known in the literature. 
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7.2 The non-planar double box 

Using the representation of ref. [31] for the two loop non-planar box (see fig. 5), we get the 
expression 



Xbox = C f 



I 



dxidX2(lX3dX4 1 5(l — X\ — X2 — X3 — Xi)dT\dT2 x 2 
(X1X3S + X 2 X 4 t c + XiX 2 Q 2 + X2X3QI ) 3+2e 



1+e 



where 



(7.21) 



Ql = (1 - Tl)(l - t 2 )s, Q 2 = tit 2 s, t c = r 2 (l - TljU+ (1 - T 2 )rit 



and 



a 



2r(3 + 2e)r(-e)r(l-e) 
r(l + e) 2 r(l-2e) 



(7.22) 



(7.23) 




Figure 5: The massless non-planar double box with all legs on-shell 

In order to avoid the singularities at the upper corners of the integration region we 
split the integral in four, mapping t\ — > 7~i/2, t\ — > 1 — t±/2 and then r 2 t 2 /2 and 
r 2 — > 1 — r 2 /2. Two of the resulting integrals can be mapped to the other two by exchanging 
xl and x$, so we end up with 



Xbox„ = C, 



43+ 



2c 



/ 



dx\dx2dx^dx4 : 5{l — x\ — x 2 — X3 — x^dT\dT2 x\^ e 

(4xiX3S + X1X2B1B2S + X2X3T1T2S + X2Xi(T2B\t + S 2 Tl«)) 3+2e 

(7.24) 



Xbox;, = C e 



43+2e 



/ 



dxidx2dx%dx4 : 5(l — X\ — x 2 — X3 — x^)dT\dT2 x\ +)L 

(4xiX 3 S + XiX2T\B 2 S + X2X3B1T2S + X 2 X 4 (i?lB 2 U + TiT 2 t)) 3+2e 

(7.25) 



where 



B 



1,2 



2 " Ti 2 



(7.26) 



Subsequently, we use the method of primary sectors on the variables xi, x 2 , X3, X4 on each 
of the above integrals, to get eight primary sectors: 

1. Xbox a i: Xbox a where x\ > x 2j 3 ) 4 has the singularity structure of x 2 + X3, and we use 
the mapping of eq. 3.12. 

2. Xbox a2 : Xbox a where x 2 > xi j3) 4 has the rather intricate singularity structure of 
the example eq. 4.45. We follow the discussion given there and decompose it to six 
integrals. 



- 33 - 



3. Xbox a 3: Xbox a where £3 > £1,2,4 has the rather intricate singularity structure of 
the example eq. 4.35. We follow the discussion given there and decompose it to five 
integrals. 

4. Xbox a 4: Xbox a where £4 > £1,2,3 has the singularity structure of eq. 4.32. We follow 
the discussion given there and decompose it in two sub-sectors, in £2, £3, each of 
which can be factorized. 

5. Xboxti: Xboxfe where x\ > £2,3,4 has the singularity structure of eq. 4.7 and we use 
the mapping of eq. 4.8. 

6. Xbox^: Xbox;, where X2 > £1,3,4 has the singularity structure of eq. 4.27 and we 
follow the mappings described there to factorize it . 

7. Xboxt3: Xboxfe where £3 > £1,2,4 has the singularity structure of eq. 4.7 and we use 
the mapping of eq. 4.8. 

8. Xbox^: Xboxfe where £4 > £1,2,3 has the singularity structure of eq. 4.1 and we use 
the mapping of eq. 4.2. 

We end up with 18 integrals to be evaluated numerically. This should be contrasted 
with the 119 sectors that are necessary if one factorizes the non-planar double box with 
sector decomposition. 

8. Conclusions 

Higher order perturbative calculations are very important for precision phenomenology at 
modern accelerator experiments. We believe that NNLO computations will be particularly 
relevant for signals of yet undiscovered physics, such as a Higgs boson or candidates of 
dark matter, in 2 — > 1 and 2 — > 2 processes. This motivates the development of powerful 
integration methods of matrix-elements of up to two virtual or real, potentially unresolved, 
partons. Such integrations entail the disentanglement of overlapping singularities. 

In this paper, we have introduced a method for the factorization of singularities based 
on non-linear transformations. As proof of principle, we presented the most singular inte- 
gral topologies which appear in NNLO double-real radiation processes with massive parti- 
cles in the final state. We find that all overlapping singularities can be factorized with our 
method, which yields a small number of numerically stable integrals. We have also applied 
our method to complicated crossed two-loop master integrals for massless QCD scatter- 
ing processes. We find that we can factorize most of the overlapping singularities with 
non-linear transformations. However, some remaining singularities are cumbersome to be 
treated purely with our method. In such situations, we employ a hybrid of our method and 
sector decomposition. This is more efficient than employing a pure sector decomposition 
approach. 

The reduction of the number of integrals which emerge in higher order corrections 
should facilitate NNLO computations. We are looking forward to applying our method for 
precision phenomenological studies of basic collider processes. 
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